######################################################################################
######################################################################################
#Replication code for:

#Mark T. Buntaine, Brigham Daniels, Colleen Devlin
#Can information outreach increase participation in community-driven development? A field experiment near Bwindi National Park, Uganda
#World Development, 106: 407-421

#Prepared by: Mark Buntaine
#Mark Buntaine contact (as of September 2018): buntaine@bren.ucsb.edu

#Compiled using R Version 3.4.3 (version "Joy in Playing") on Mac running OS X 10.13.6
#+ fig.width=6.5, fig.height=8.25
######################################################################################
######################################################################################

root <- "~/Google Drive/Bwindi project/Code/WD_Replication"
#Note: all code and data files should be placed in a single directory
#Note: change the file path "root" to the local working directory
setwd(root)

######################################################################################
###Packages
######################################################################################
library(interplot) #Version 0.1.5
library(interflex) #Version 1.0.3
library(lfe) #Version 2.6-2291
library(Lmoments) #Version 1.2-3, for inter.binning.90 function


######################################################################################
###Functions
######################################################################################

mean.imp <- function(data,var,cluster){
  for (i in 1:length(unique(data[,cluster]))){
    imp.value <- mean(data[data[,cluster]==unique(data[,cluster])[i], var], na.rm=T)
    data[data[,cluster]==unique(data[,cluster])[i] & is.na(data[,var]), var] <- imp.value
  }
  return(data)
} #This function automates mean imputation by cluster


felm.ri <- function(formula, dta, treat.var, sims, ...){
  require(lfe)
  
  ate <- coef(felm(formula, data=dta))[1]
  N <- felm(formula, data=dta)$N
  ate.samp.dist <- rep(NA,sims)
  
  for (i in 1:sims){
    dta[,treat.var] <- dta[,ncol(M)+i]
    ate.samp.dist[i] <- coef(felm(formula, data=dta))[1]
  }
  
  p.two.way <- sum(abs(ate)<abs(ate.samp.dist))/sims
  p.one.way.greater <- sum(ate<ate.samp.dist)/sims
  p.one.way.lesser <- sum(ate>ate.samp.dist)/sims
  se <- sd(ate.samp.dist)
  
  coun <- list("ate" = ate, "ate.samp.dist" = ate.samp.dist, "se"=se, "p.two.way" = p.two.way, "p.one.way.greater" = p.one.way.greater, "p.one.way.lesser" = p.one.way.lesser, "N" = N)
  return(coun)
}


lm.ri <- function(formula, dta, treat.var, sims, ...){
  require(lfe)
  
  ate <- coef(lm(formula, data=dta))[2]
  N <- nobs(lm(formula, data=dta))
  ate.samp.dist <- rep(NA,sims)
  
  for (i in 1:sims){
    dta[,treat.var] <- dta[,ncol(M)+i]
    ate.samp.dist[i] <- coef(lm(formula, data=dta))[2]
  }
  
  p.two.way <- sum(abs(ate)<abs(ate.samp.dist))/sims
  p.one.way.greater <- sum(ate<ate.samp.dist)/sims
  p.one.way.lesser <- sum(ate>ate.samp.dist)/sims
  se <- sd(ate.samp.dist)
  
  coun <- list("ate" = ate, "ate.samp.dist" = ate.samp.dist, "se"=se, "p.two.way" = p.two.way, "p.one.way.greater" = p.one.way.greater, "p.one.way.lesser" = p.one.way.lesser, "N" = N)
  return(coun)
}

source('inter.binning.90.R')
source('multiplot.R')

#Function to reset par() show that blocks can be run sequentially and without a plotting device
default.par <- par()


######################################################################################
###Data input
######################################################################################

B <- read.csv("baseline.csv", stringsAsFactors=FALSE)
E <- read.csv("endline.csv", stringsAsFactors=FALSE)
endline.subjects <- read.csv("endline_call_list.csv", row.names=1, stringsAsFactors=FALSE)[,1:4]
recruited.subjects <- read.csv("recruited_subjects.csv", stringsAsFactors=FALSE)
smsid <- read.csv("sms_received_list.csv", stringsAsFactors=FALSE)
Responsiveness <- read.csv("digital_responses.csv", stringsAsFactors=FALSE)
con.mat <- read.csv("contiguity_mat1.csv", stringsAsFactors=FALSE)
con2.mat <- read.csv("contiguity_mat2.csv", stringsAsFactors=FALSE)
con3.mat <- read.csv("contiguity_mat3.csv", stringsAsFactors=FALSE)
#Note the "con.mat" files have village "Kiziiba" manually removed because merged with "Buguro" at cleaning phase
spill.dist <- read.csv("spillover_matrix.csv", row.names=1)
loc.master <- read.csv("location_master.csv", stringsAsFactors=FALSE)
lc1.planning <- read.csv("planning_meeting_attendance.csv", stringsAsFactors=FALSE)
sbj_RI <- read.csv("assignment_ri.csv", stringsAsFactors=FALSE)


######################################################################################
###Data setup
######################################################################################

#Removing duplicate subject.id's in endline.subjects
endline.dup <- endline.subjects$Subject_ID_number[duplicated(endline.subjects$Subject_ID_number)]
see <- endline.subjects[endline.subjects$Subject_ID_number %in% endline.dup,]
see <- see[order(see$Subject_ID_number),] #Visual inspection: all subject id's are exact duplicates
M <- endline.subjects[!duplicated(endline.subjects$Subject_ID_number),] #Removing duplicated subject id's
length(unique(M$Subject_ID_number)) #This is the merge point, as all subject id's are unique

#Preparing baseline merge
B <- subset(B, B_subject_refuse_to_participate=="no") #Removing subjects who refused, no data
length(unique(B$Subject_ID_number)) #There are several duplicate subject id's in the raw baseline file
B.duplicated.list <- B$Subject_ID_number[duplicated(B$Subject_ID_number)] #List of duplicate subject id's in raw endline file
see <- B[B$Subject_ID_number %in% B.duplicated.list,]
see <- see[order(see$Subject_ID_number),]
B <- B[!(B$Subject_ID_number %in% B.duplicated.list),]
M <- merge(M,B, by="Subject_ID_number", all.x=TRUE) #still 1723 observations

#Preparing endline merge
E <- subset(E, E_1_Did_the_subject_refuse_to_p=="no") #Removing subjects who refused, no data
E.duplicated.list <- E$Subject_ID_number[duplicated(E$Subject_ID_number)] #List of duplicate subject id's in raw endline file
see <- E[E$Subject_ID_number %in% E.duplicated.list,]
see <- see[order(see$Subject_ID_number),]
E <- E[!(E$Subject_ID_number %in% E.duplicated.list),]
M <- merge(M,E, by="Subject_ID_number", all.x=TRUE) #still 1723 observations

#Merging treatment status and block from the assignment file
treat.dta <- recruited.subjects[,c(1,10,11,12)] #Getting subcounty, treatment status, and block
M <- merge(M,treat.dta, by="Subject_ID_number", all.x=TRUE) #still 1723 observations
M$subcounty.blocking <- ifelse(M$subcounty.blocking=="Rutenga", "Singles", M$subcounty.blocking) #Because of exclusion, Rutenga subcounty has only one village, moving to "Singles" block

#TextIt data
m7<-table(smsid$Subject_ID_number)
m7<-data.frame(m7)
names(m7)<-c("Subject_ID_number","participate.m7")
M <- merge(M,m7, by="Subject_ID_number", all.x=TRUE) #still 1723 observations

#Location cleaning
villages <- loc.master$village.updated
for (i in 1:length(villages)){
  villages[i] <- ifelse(loc.master$village.cleaned[i]=="", loc.master$village.updated[i], loc.master$village.cleaned[i])
}
villages[duplicated(villages)] #Checking for errors; looks good though some original errors resulted in villages with mixed treatment assignment (Kajumo, Omumbuga)

see <- subset(loc.master,excluded==0)
for (i in 1:length(villages)){
  villages[i] <- ifelse(see$village.cleaned[i]=="", see$village.updated[i], see$village.cleaned[i])
}
villages[duplicated(villages)] #Checking for errors; looks good though some original errors resulted in villages with mixed treatment assignment (Kajumo, Omumbuga)
see$final.id <- ifelse(see$location.id.cleaned=="",see$location.id,see$location.id.cleaned)
unique(see$final.id)

#unique(M$Updated.Village)[!(unique(M$Updated.Village) %in% villages)]
M$Updated.Village <- ifelse(M$Updated.Village=="Katooma","Katoma",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Kijumwa","Kijuma",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Kyogo (Kabale)","Kyogo (Rubanda)",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Mushija","Musheija",M$Updated.Village)
M$Updated.Village <- ifelse(M$Updated.Village=="Ndeego","Ndego",M$Updated.Village)

#Adding variables to probe spillover
spill.village.list <- names(spill.dist)
spill.village.list[!(spill.village.list %in% recruited.subjects$updated.village)] #list in spillover matrix that don't match

#con.village.list <- con.mat$Village..m.
#con.village.list[!(con.village.list %in% M$Updated.Village)] #list in contiguity matrix that don't match M village names
#con.village.list[!(con.village.list %in% loc.master$village.updated)] #list in contiguity matrix that don't match location master original names

#Cleaning the con.mat file for village names, see email thread "village names in matrices" from Jan-2017
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kashija/ Rubuguri","Kashija",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Higabiro/ Rubuguri","Higabiro (Rubuguri)",con.mat$Village..m.)
#Higabiro (Rubuguri): this is an RS village, but we do not have any subjects from this village
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nteeko","Nteko",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakabungo( Rubimbwa)","Nyakabungo (Rubimbwa)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakabungo(Bushura)","Nyakabungo (Bushura)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Ishaya (Mucusye)","Ishaya",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kiziba","Kiziiba",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Bushegyenyi","Bushegenyi",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kyibingo","Kibingo",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kitahulira","Kitahurira (Hakikome)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakabungo","Nyakabungo (Kashekyera)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Buzaniza","Buzaniro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Mburamaizi","Mburameizi",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Bishanyu","Bishayu",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Katooma","Katoma",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Rwensaniro","Rwensanziro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Rugando","Rugandu",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Ntunagamo","Ntungamo",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Nyakaranaga","Nyakaranga",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kacwamuhoro","Kachwamuhoro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kitahurira","Kitahurira (Kashasha)",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Ndeego","Ndego",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="kagogo","Kagogo",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Byamihanda","Ryamihanda",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Mushija","Musheija",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kijumwa","Kijuma",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kiziiba","Bugoro",con.mat$Village..m.)
con.mat$Village..m. <- ifelse(con.mat$Village..m.=="Kyogo (Kabale)","Kyogo (Rubanda)",con.mat$Village..m.)

#Cleaning the con2.mat file for village names, see email thread "village names in matrices" from Jan-2017
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kashija/ Rubuguri","Kashija",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Higabiro/ Rubuguri","Higabiro (Rubuguri)",con2.mat$Village..m.)
#Higabiro (Rubuguri): this is an RS village, but we do not have any subjects from this village
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nteeko","Nteko",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakabungo( Rubimbwa)","Nyakabungo (Rubimbwa)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakabungo(Bushura)","Nyakabungo (Bushura)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Ishaya (Mucusye)","Ishaya",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kiziba","Kiziiba",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Bushegyenyi","Bushegenyi",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kyibingo","Kibingo",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kitahulira","Kitahurira (Hakikome)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakabungo","Nyakabungo (Kashekyera)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Buzaniza","Buzaniro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Mburamaizi","Mburameizi",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Bishanyu","Bishayu",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Katooma","Katoma",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Rwensaniro","Rwensanziro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Rugando","Rugandu",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Ntunagamo","Ntungamo",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Nyakaranaga","Nyakaranga",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kacwamuhoro","Kachwamuhoro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kitahurira","Kitahurira (Kashasha)",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Ndeego","Ndego",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="kagogo","Kagogo",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Byamihanda","Ryamihanda",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Mushija","Musheija",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kijumwa","Kijuma",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kiziiba","Bugoro",con2.mat$Village..m.)
con2.mat$Village..m. <- ifelse(con2.mat$Village..m.=="Kyogo (Kabale)","Kyogo (Rubanda)",con2.mat$Village..m.)

#Cleaning the con3.mat file for village names, see email thread "village names in matrices" from Jan-2017
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kashija/ Rubuguri","Kashija",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Higabiro/ Rubuguri","Higabiro (Rubuguri)",con3.mat$Village..m.)
#Higabiro (Rubuguri): this is an RS village, but we do not have any subjects from this village
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nteeko","Nteko",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakabungo( Rubimbwa)","Nyakabungo (Rubimbwa)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakabungo(Bushura)","Nyakabungo (Bushura)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Ishaya (Mucusye)","Ishaya",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kiziba","Kiziiba",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Bushegyenyi","Bushegenyi",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kyibingo","Kibingo",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kitahulira","Kitahurira (Hakikome)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakabungo","Nyakabungo (Kashekyera)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Buzaniza","Buzaniro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Mburamaizi","Mburameizi",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Bishanyu","Bishayu",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Katooma","Katoma",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Rwensaniro","Rwensanziro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Rugando","Rugandu",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Ntunagamo","Ntungamo",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Nyakaranaga","Nyakaranga",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kacwamuhoro","Kachwamuhoro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kitahurira","Kitahurira (Kashasha)",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Ndeego","Ndego",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="kagogo","Kagogo",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Byamihanda","Ryamihanda",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Mushija","Musheija",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kijumwa","Kijuma",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kiziiba","Bugoro",con3.mat$Village..m.)
con3.mat$Village..m. <- ifelse(con3.mat$Village..m.=="Kyogo (Kabale)","Kyogo (Rubanda)",con3.mat$Village..m.)

#Cleaning loc.master names into one vector
loc.master$village.final <- ifelse(loc.master$village.cleaned=="", loc.master$village.updated, loc.master$village.cleaned)
loc.master$village.final <- ifelse(loc.master$village.final=="Kikangaga ", "Kikangaga", loc.master$village.final)

#Checking the location.id cleaning in all files
con.mat$Village..m.[!(con.mat$Village..m. %in% loc.master$village.final)] #Notes: only Higabiro (Rubuguri), no subjects in study
loc.master$village.final[!(loc.master$village.final %in% con.mat$Village..m.)] #Notes: "Kinaba" "Kinaba" "Kyogo"  "Rukere" all excluded
unique(M$Updated.Village) %in% loc.master$village.final #Note: all appear in cleaned loc.master
unique(M$Updated.Village)[!(unique(M$Updated.Village) %in% con.mat$Village..m.)] #"Kinaba" does not specify which one, remove Kinaba from dataframe?

#Correcting Kinaba village names (redistricted), for those we were able to reach, see "Corrected villages.xls" and email thread "Bwindi Village Confusion"
M$Updated.Village[M$Subject_ID_number %in% c(2581,2583,2587)] <- "Nyabiha"
M$Updated.Village[M$Subject_ID_number %in% c(2518)] <- "Kyabworo"
M <- M[M$Updated.Village!="Kinaba",] #Removing "Kinaba" from dataframe per exclusion from revenue sharing during '16/'17 round

#Creating spillover indicator in M
treated.villages <- loc.master$village.final[loc.master$treat==1]
#Note: rows/columns in con.mat must be ordered identically, code below checks
#col.n <- colnames(con.mat)[-1]
#ex <- dataframe(con.mat$Village..m.,col.n)
#con.mat$Village..m.[48] <- "Kiziba" #Correction discussed in 2/20/17 email to Colleen

for (i in 1:nrow(M)){
  row.test <- con.mat[M$Updated.Village[i]==con.mat$Village..m.,-1]==1
  row.test2 <- con2.mat[M$Updated.Village[i]==con2.mat$Village..m.,-1]==1
  row.test3 <- con3.mat[M$Updated.Village[i]==con3.mat$Village..m.,-1]==1
  con.villages <- con.mat$Village..m.[row.test]
  con2.villages <- con2.mat$Village..m.[row.test2]
  con3.villages <- con3.mat$Village..m.[row.test3]
  M$S1[i] <- sum(con.villages %in% treated.villages)
  M$S2[i] <- sum(con2.villages %in% treated.villages)
  M$S3[i] <- sum(con3.villages %in% treated.villages)
  
}
M$S1 <- as.factor(M$S1)
M$S2 <- as.factor(M$S2)
M$S3 <- as.factor(M$S3)

#Treatment Density (number of subjects per village)
for (i in 1:nrow(M)){
  sub <- subset(M, Updated.Village==M$Updated.Village[i])
  M$subjects.per.village[i] <- nrow(sub)
}

#Attrition
M$attr <- ifelse(is.na(M$E_id), 1, 0)

#Setting up RI object
source('outcome_var_mean_imp_setup.R')
sbj_RI <- sbj_RI[,c(1,5:ncol(sbj_RI))] #selecting only the treatment status to avoid introducing duplicate data to M


######################################################################################
###Table 1: Demographic data collected during endline survey by treatment group
######################################################################################
#Note: the table is formatted manually based on the values computed below

#Gender
prop.table(table(M$treat,M$E_subject_s_gender),1)

#Age
summary(M$E_Subject_s_age[M$treat==1])
summary(M$E_Subject_s_age[M$treat==0])

#Partial or full literacy
table(M$E_Can_you_read_001[M$treat==1])
table(M$E_Can_you_read_001[M$treat==0])

#Monthly income
table(M$E_Subject_s_income[M$treat==1])
table(M$E_Subject_s_income[M$treat==0])

#Subjects
table(M$treat)

#Endline Survey Participants (subjects in column without attrition)
table(M$treat,M$attr)


######################################################################################
###Figure 3: Unadjusted means for survey-based outcomes measures by treatment condition
######################################################################################

setwd("~/Google Drive/Bwindi project/Code/WD_Replication")
source('outcome_var_mean_imp_setup.R')

set.seed(101)
sims <- 10000

##Support for Conservation of Bwindi
b.c.b6.t <- mean(M$b.conservation.b6[M$treat==1 & !is.na(M$e.conservation.b6)], na.rm=T)
b.c.b6.c <- mean(M$b.conservation.b6[M$treat==0 & !is.na(M$e.conservation.b6)], na.rm=T)

c.b6.t <- mean(M$e.conservation.b6[M$treat==1], na.rm=T)
c.b6.c <- mean(M$e.conservation.b6[M$treat==0], na.rm=T)

c.b6.t.sims <- rep(NA, sims)
c.b6.c.sims <- rep(NA, sims)

for (i in 1:sims){
  c.b6.t.sims[i] <- mean(sample(M$e.conservation.b6[M$treat==1], replace=T), na.rm=T)
  c.b6.c.sims[i] <- mean(sample(M$e.conservation.b6[M$treat==0], replace=T), na.rm=T)
}

##Satisfaction with Revenue Sharing
b.s.b3.t <- mean(M$b.satisfaction.b3[M$treat==1  & !is.na(M$e.satisfaction.b3)], na.rm=T)
b.s.b3.c <- mean(M$b.satisfaction.b3[M$treat==0 & !is.na(M$e.satisfaction.b3)], na.rm=T)

s.b3.t <- mean(M$e.satisfaction.b3[M$treat==1], na.rm=T)
s.b3.c <- mean(M$e.satisfaction.b3[M$treat==0], na.rm=T)

s.b3.t.sims <- rep(NA, sims)
s.b3.c.sims <- rep(NA, sims)

for (i in 1:sims){
  s.b3.t.sims[i] <- mean(sample(M$e.satisfaction.b3[M$treat==1], replace=T), na.rm=T)
  s.b3.c.sims[i] <- mean(sample(M$e.satisfaction.b3[M$treat==0], replace=T), na.rm=T)
}

##Satisfaction with Park Management
b.s.b2.t <- mean(M$b.satisfaction.b2[M$treat==1 & !is.na(M$e.satisfaction.b2)], na.rm=T)
b.s.b2.c <- mean(M$b.satisfaction.b2[M$treat==0 & !is.na(M$e.satisfaction.b2)], na.rm=T)

s.b2.t <- mean(M$e.satisfaction.b2[M$treat==1], na.rm=T)
s.b2.c <- mean(M$e.satisfaction.b2[M$treat==0], na.rm=T)

s.b2.t.sims <- rep(NA, sims)
s.b2.c.sims <- rep(NA, sims)

for (i in 1:sims){
  s.b2.t.sims[i] <- mean(sample(M$e.satisfaction.b2[M$treat==1], replace=T), na.rm=T)
  s.b2.c.sims[i] <- mean(sample(M$e.satisfaction.b2[M$treat==0], replace=T), na.rm=T)
}

##Know Right Person to Contact about RS
b.s.b12.t <- mean(M$b.participate.b12[M$treat==1 & !is.na(M$e.participate.b12)], na.rm=T)
b.s.b12.c <- mean(M$b.participate.b12[M$treat==0 & !is.na(M$e.participate.b12)], na.rm=T)

s.b12.t <- mean(M$e.participate.b12[M$treat==1], na.rm=T)
s.b12.c <- mean(M$e.participate.b12[M$treat==0], na.rm=T)

s.b12.t.sims <- rep(NA, sims)
s.b12.c.sims <- rep(NA, sims)

for (i in 1:sims){
  s.b12.t.sims[i] <- mean(sample(M$e.participate.b12[M$treat==1], replace=T), na.rm=T)
  s.b12.c.sims[i] <- mean(sample(M$e.participate.b12[M$treat==0], replace=T), na.rm=T)
}

##Opportunities to Communicate with UWA about RS
b.s.b11.t <- mean(M$b.participate.b11[M$treat==1 & !is.na(M$e.participate.b11)], na.rm=T)
b.s.b11.c <- mean(M$b.participate.b11[M$treat==0 & !is.na(M$e.participate.b11)], na.rm=T)

s.b11.t <- mean(M$e.participate.b11[M$treat==1], na.rm=T)
s.b11.c <- mean(M$e.participate.b11[M$treat==0], na.rm=T)

s.b11.t.sims <- rep(NA, sims)
s.b11.c.sims <- rep(NA, sims)

for (i in 1:sims){
  s.b11.t.sims[i] <- mean(sample(M$e.participate.b11[M$treat==1], replace=T), na.rm=T)
  s.b11.c.sims[i] <- mean(sample(M$e.participate.b11[M$treat==0], replace=T), na.rm=T)
}

##Participated in RS meeting in past month
M$e.participate.e4 <- as.numeric(M$e.participate.e4)
e4.t <- mean(M$e.participate.e4[M$treat==1], na.rm=T)
e4.c <- mean(M$e.participate.e4[M$treat==0], na.rm=T)

e4.t.sims <- rep(NA, sims)
e4.c.sims <- rep(NA, sims)

for (i in 1:sims){
  e4.t.sims[i] <- mean(sample(M$e.participate.e4[M$treat==1], replace=T), na.rm=T)
  e4.c.sims[i] <- mean(sample(M$e.participate.e4[M$treat==0], replace=T), na.rm=T)
}

##Opportunities to Participate in RS
b.s.b10.t <- mean(M$b.participate.b10[M$treat==1 & !is.na(M$e.participate.b10)], na.rm=T)
b.s.b10.c <- mean(M$b.participate.b10[M$treat==0 & !is.na(M$e.participate.b10)], na.rm=T)

s.b10.t <- mean(M$e.participate.b10[M$treat==1], na.rm=T)
s.b10.c <- mean(M$e.participate.b10[M$treat==0], na.rm=T)

#t.test(x=M$e.participate.b10[M$treat==1], y=M$e.participate.b10[M$treat==0], na.action=na.omit)

s.b10.t.sims <- rep(NA, sims)
s.b10.c.sims <- rep(NA, sims)

for (i in 1:sims){
  s.b10.t.sims[i] <- mean(sample(M$e.participate.b10[M$treat==1], replace=T), na.rm=T)
  s.b10.c.sims[i] <- mean(sample(M$e.participate.b10[M$treat==0], replace=T), na.rm=T)
}

##Know How RS Works in my Village
b.b9.t <- mean(M$b.knowledge.b9[M$treat==1 & !is.na(M$e.knowledge.b9)], na.rm=T)
b.b9.c <- mean(M$b.knowledge.b9[M$treat==0 & !is.na(M$e.knowledge.b9)], na.rm=T)

b9.t <- mean(M$e.knowledge.b9[M$treat==1], na.rm=T)
b9.c <- mean(M$e.knowledge.b9[M$treat==0], na.rm=T)

#t.test(x=M$e.knowledge.b9[M$treat==1], y=M$e.knowledge.b9[M$treat==0], na.action=na.omit)

b9.t.sims <- rep(NA, sims)
b9.c.sims <- rep(NA, sims)

for (i in 1:sims){
  b9.t.sims[i] <- mean(sample(M$e.knowledge.b9[M$treat==1], replace=T), na.rm=T)
  b9.c.sims[i] <- mean(sample(M$e.knowledge.b9[M$treat==0], replace=T), na.rm=T)
}

##Know How RS Works Generally
b.b8.t <- mean(M$b.knowledge.b8[M$treat==1 & !is.na(M$e.knowledge.b8)], na.rm=T)
b.b8.c <- mean(M$b.knowledge.b8[M$treat==0 & !is.na(M$e.knowledge.b8)], na.rm=T)

b8.t <- mean(M$e.knowledge.b8[M$treat==1], na.rm=T)
b8.c <- mean(M$e.knowledge.b8[M$treat==0], na.rm=T)

b8.t.sims <- rep(NA, sims)
b8.c.sims <- rep(NA, sims)

for (i in 1:sims){
  b8.t.sims[i] <- mean(sample(M$e.knowledge.b8[M$treat==1], replace=T), na.rm=T)
  b8.c.sims[i] <- mean(sample(M$e.knowledge.b8[M$treat==0], replace=T), na.rm=T)
}

##Know How RS Works in my Village
b.b7.t <- mean(M$b.knowledge.b7[M$treat==1 & !is.na(M$e.knowledge.b7)], na.rm=T)
b.b7.c <- mean(M$b.knowledge.b7[M$treat==0 & !is.na(M$e.knowledge.b7)], na.rm=T)

b7.t <- mean(M$e.knowledge.b7[M$treat==1], na.rm=T)
b7.c <- mean(M$e.knowledge.b7[M$treat==0], na.rm=T)

b7.t.sims <- rep(NA, sims)
b7.c.sims <- rep(NA, sims)

for (i in 1:sims){
  b7.t.sims[i] <- mean(sample(M$e.knowledge.b7[M$treat==1], replace=T), na.rm=T)
  b7.c.sims[i] <- mean(sample(M$e.knowledge.b7[M$treat==0], replace=T), na.rm=T)
}

#Plots (with baseline values, ordered and color-coded)
#png("Bwindi_ParticipationSurveyMeans_170502.png",
#    width=6.5,height=8.25,units="in", res=300)
par(oma=c(0.5,0.5,0.5,0.5),mar=c(3.4,2.3,0.5,1.9), mgp=c(2.2,0.8,0), bty="n", cex.lab=1.2, cex.axis=1.2)
layout(matrix(1:11, ncol = 1), widths = 1, heights = c(rep(1,10),0.4), respect = FALSE)

plot(y=c(2.5,1.5),x=c(b8.t,b8.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,4), yaxt="n", ylab="", xlab="Know How RS Works Generally", col.lab="darkgreen", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4), labels=c("Disagree","Somewhat Disagree","Somewhat Agree","Agree"))
lines(y=c(2.5,2.5),x=c(sort(b8.t.sims)[250],sort(b8.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(b8.c.sims)[250],sort(b8.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.b8.t,b.b8.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(b9.t,b9.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,4), yaxt="n", ylab="", xlab="Know How RS Works in my Village", col.lab="darkgreen", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4), labels=c("Disagree","Somewhat Disagree","Somewhat Agree","Agree"))
lines(y=c(2.5,2.5),x=c(sort(b9.t.sims)[250],sort(b9.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(b9.c.sims)[250],sort(b9.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.b9.t,b.b9.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(s.b12.t,s.b12.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,4), yaxt="n", ylab="", xlab="Know Right Person to Contact about RS", col.lab="darkgreen", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4), labels=c("Disagree","Somewhat Disagree","Somewhat Agree","Agree"))
lines(y=c(2.5,2.5),x=c(sort(s.b12.t.sims)[250],sort(s.b12.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(s.b12.c.sims)[250],sort(s.b12.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.s.b12.t,b.s.b12.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(s.b11.t,s.b11.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,5), yaxt="n", ylab="", xlab="Opportunities to Communicate with UWA about RS", col.lab="red", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4,5), labels=c("Very Dissatisfied","Somewhat Dissatisfied","Neutral","Somewhat Satisfied","Very Satisfied"))
lines(y=c(2.5,2.5),x=c(sort(s.b11.t.sims)[250],sort(s.b11.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(s.b11.c.sims)[250],sort(s.b11.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.s.b11.t,b.s.b11.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(s.b10.t,s.b10.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,4), yaxt="n", ylab="", xlab="Opportunities to Participate in Revenue Sharing", col.lab="red", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4), labels=c("Disagree","Somewhat Disagree","Somewhat Agree","Agree"))
lines(y=c(2.5,2.5),x=c(sort(s.b10.t.sims)[250],sort(s.b10.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(s.b10.c.sims)[250],sort(s.b10.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.s.b10.t,b.s.b10.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(e4.t,e4.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(0,1), yaxt="n", ylab="", xlab="Participated in RS meeting", col.lab="red", xaxt="n", font.lab=2)
axis(1, at=c(0,1), labels=c("No","Yes"))
lines(y=c(2.5,2.5),x=c(sort(e4.t.sims)[250],sort(e4.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(e4.c.sims)[250],sort(e4.c.sims)[9750]), lwd=1.5)

plot(y=c(2.5,1.5),x=c(b7.t,b7.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,5), yaxt="n", ylab="", xlab="Satisfied with Information from UWA about RS", col.lab="blue", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4,5), labels=c("Very Dissatisfied","Somewhat Dissatisfied","Neutral","Somewhat Satisfied","Very Satisfied"))
lines(y=c(2.5,2.5),x=c(sort(b7.t.sims)[250],sort(b7.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(b7.c.sims)[250],sort(b7.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.b7.t,b.b7.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(s.b3.t,s.b3.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,5), yaxt="n", ylab="", xlab="Satisfied with Revenue Sharing", col.lab="blue", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4,5), labels=c("Very Dissatisfied","Somewhat Dissatisfied","Neutral","Somewhat Satisfied","Very Satisfied"))
lines(y=c(2.5,2.5),x=c(sort(s.b3.t.sims)[250],sort(s.b3.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(s.b3.c.sims)[250],sort(s.b3.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.s.b3.t,b.s.b3.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(s.b2.t,s.b2.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5),xlim=c(1,5), yaxt="n", ylab="", xlab="Satisfied with Park Management", col.lab="blue", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4,5), labels=c("Very Dissatisfied","Somewhat Dissatisfied","Neutral","Somewhat Satisfied","Very Satisfied"))
lines(y=c(2.5,2.5),x=c(sort(s.b2.t.sims)[250],sort(s.b2.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(s.b2.c.sims)[250],sort(s.b2.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.s.b2.t,b.s.b2.c), pch=c(16,1), cex=1, col="red")

plot(y=c(2.5,1.5),x=c(c.b6.t,c.b6.c), pch=c(16,1), cex=2, ylim=c(0.5,3.5), xlim=c(1,4), yaxt="n", ylab="", xlab="Importance of Protecting Bwindi", col.lab="blue", xaxt="n", font.lab=2)
axis(1, at=c(1,2,3,4), labels=c("Not At All Important","Not Very Important","Somewhat Important","Very Important"))
lines(y=c(2.5,2.5),x=c(sort(c.b6.t.sims)[250],sort(c.b6.t.sims)[9750]), lwd=1.5)
lines(y=c(1.5,1.5),x=c(sort(c.b6.c.sims)[250],sort(c.b6.c.sims)[9750]), lwd=1.5)
points(y=c(2.5,1.5), x=c(b.c.b6.t,b.c.b6.c), pch=c(16,1), cex=1, col="red")

par(mar=c(0,0.1,0.1,0.1))
plot(1, type="n", xlim=c(0,10), ylim=c(0,2), xaxt="n", yaxt="n", ylab="", xlab="")
rect(xleft=1.7,xright=8.3,ybottom=0,ytop=2, col="azure2")
points(x=c(2.3,6,2.3,6),y=c(1.45,1.45,0.35,0.35), pch=c(16,21,16,21), col=c("black","black","red","red"), bg="white", cex=c(2,2,1,1))
text(x=2.4, y=1.35, pos=4, "Treatment (Endline)", cex=1.2)
text(x=6.1, y=1.35, pos=4, "Control (Endline)", cex=1.2)
text(x=2.4, y=0.35, pos=4, "Treatment (Baseline)", cex=1)
text(x=6.1, y=0.35, pos=4, "Control (Baseline)", cex=1)

#dev.off()


######################################################################################
###Figure 4: Treatment effects on survey items
######################################################################################

setwd("~/Google Drive/Bwindi project/Code/WD_Replication")
source('outcome_var_mean_imp_setup.R')
M.ri <- merge(M, sbj_RI, by="Subject_ID_number", all.x=TRUE)


#Model-assisted RI
mod.h7.ri <- felm.ri(e.knowledge.b7 ~ treat + b.knowledge.b7 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                     dta=M.ri, treat.var = "treat", sims=10000)

mod.h7.2.ri <- felm.ri(e.knowledge.b8 ~ treat + b.knowledge.b8 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

mod.h7.3.ri <- felm.ri(e.knowledge.b9 ~ treat + b.knowledge.b9 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

mod.h8.1.ri <- felm.ri(e.participate.b10 ~ treat + b.participate.b10 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

mod.h8.3.ri <- felm.ri(e.participate.e4 ~ treat + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

mod.h9.1.ri <- felm.ri(e.participate.b11 ~ treat + b.participate.b11 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                       dta=M.ri, treat.var = "treat", sims=10000)

mod.h9.2.ri <- lm.ri(e.participate.b12 ~ treat + b.participate.b12 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M.ri, treat.var = "treat", sims=10000)

nod.h10mgmt.ri <- felm.ri(e.satisfaction.b2 ~ treat + b.satisfaction.b2 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                          dta=M.ri, treat.var = "treat", sims=10000)

nod.h11.ri <- felm.ri(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village | subcounty.blocking | 0 | subcounty.blocking,
                      dta=M.ri, treat.var = "treat", sims=10000)

nod.h13.ri <- lm.ri(e.conservation.b6 ~ treat + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                    dta=M.ri, treat.var = "treat", sims=10000)

#Plot
names <- c("Importance of Protecting Bwindi","Satisfied with Park Management","Satisfied with Revenue Sharing","Satisfied with Information from UWA about RS","Participated in RS meeting (past month)","Opportunities to Participate in RS","Opportunities to Communicate with UWA about RS","Know Right Person to Contact about RS","Know How RS Works in my Village","Know How RS Works Generally")
ates <- c(nod.h13.ri$ate,nod.h10mgmt.ri$ate,nod.h11.ri$ate,mod.h7.ri$ate,mod.h8.3.ri$ate,mod.h8.1.ri$ate,mod.h9.1.ri$ate,mod.h9.2.ri$ate,mod.h7.3.ri$ate,mod.h7.2.ri$ate)
ses <- c(nod.h13.ri$se,nod.h10mgmt.ri$se,nod.h11.ri$se,mod.h7.ri$se,mod.h8.3.ri$se,mod.h8.1.ri$se,mod.h9.1.ri$se,mod.h9.2.ri$se,mod.h7.3.ri$se,mod.h7.2.ri$se)
dta.coef <- data.frame(names,ates,ses)
dta.coef$ylo <- dta.coef$ates - 1.645*dta.coef$ses
dta.coef$yhi <- dta.coef$ates + 1.645*dta.coef$ses
dta.coef$names <- factor(dta.coef$names, levels=dta.coef$names)
dta.coef$clr <- c(rep("blue",4),rep("red",3),rep("darkgreen",3))

p2 = ggplot(dta.coef, aes(x=names, y=ates, ymin=ylo, ymax=yhi)) + geom_errorbar(colour=dta.coef$clr, width=0.2, size=1.4) + 
  theme_bw(base_size=16) + coord_flip() + geom_hline(aes(yintercept=0), lty=2) + ylab("Treatment Effects on Survey Items") + xlab("") + geom_point(size=4, colour=dta.coef$clr) +
  theme(axis.title = element_text(size=14))
p2 #Copied at 800x420


######################################################################################
###Figure 5: Digital engagement
######################################################################################

qs<-Responsiveness[, c("H_6.3", "H_6.10", "H_6.17", "H_6.24","H_7.1", "H_7.8","H_7.13", "H_8.25", "H_8.26", "H_9.9", "H_9.16", "H_9.23","H_9.30")]
row.names(qs)<-c("Treat", "Control")

names(qs)<-c("Q1", "Q2", "Q3", "Q4", "Q5", "Q6", "Q7", "Q8", "Q9", "Q10", "Q11", "Q12", "Q13")  
qs.names <- names(qs)

par(default.par)
par(mar = c(4.1, 4.1, 2.1, 2.1))
treat.prop <- qs[1,] / table(M$treat)[2]
control.prop <- qs[2,] / table(M$treat)[1]
x <- 1:13
plot(x=x, y=treat.prop, type="l", col="red", lwd=3, xlab="Question Number", ylab="Proportion Responding", ylim=c(0,0.2), xaxt="n")
axis(1, at=1:13, labels=qs.names)
lines(x=x, y=control.prop, col="grey", lwd=3)
legend("topright",legend=c("Placebo","Treated"), col=c("grey","red"), fill=c("grey","red"))


######################################################################################
###Figure 6: Digital engagement
######################################################################################

recruited.subjects$location.id <- paste(recruited.subjects$Subcounty,recruited.subjects$updated.village,sep=".")
lc1.planning$location.id <- paste(lc1.planning$Sub.County,lc1.planning$Village,sep=".")

loc.check <- lc1.planning$location.id
loc.check %in% recruited.subjects$location.id
see <- lc1.planning[!(loc.check %in% recruited.subjects$location.id),] #LC1 data without matches
#Kirundo.Higabiro: no apparent match with sampling frame
#Kayonza.Musheija
lc1.planning$location.id <- ifelse(lc1.planning$location.id=="Kayonza.Musheija", "Kayonza.Mushija", lc1.planning$location.id)
#Mpungu.Hakikome (Kitahurira): Parish/Village switched, see below
#Kirima.Mutojo/Bugandaro: combined meeting for unexplained reasons, excluding from analysis
#Ruhija.Kyogo (Rubanda): no apparent match with sampling frame

unique(recruited.subjects$location.id) %in% loc.check
unique(recruited.subjects$location.id)[!(unique(recruited.subjects$location.id) %in% loc.check)]

recruited.subjects$location.id <- ifelse(recruited.subjects$location.id=="Mpungu.Kitahurira (Hakikome)", "Mpungu.Hakikome (Kitahurira)", recruited.subjects$location.id)

for (i in 1:nrow(lc1.planning)){
  sub <- subset(recruited.subjects, location.id==lc1.planning$location.id[i])
  lc1.planning$treat[i] <- sub$treat[1]
}

lc1.planning$Reason.Attendance.Data.Missing[c(5,19,25,35,47,48,55,75,76,81,85,94)] <- "Subcounty project"
#Note: this is cleaning the column

lc1.planning.final <- subset(lc1.planning, !is.na(treat))
sum(!is.na(lc1.planning.final$Attendance))

lc1.planning.final$attr <- ifelse(is.na(lc1.planning.final$Attendance), 1, 0)

attr.mod <- lm(attr ~ treat, data=lc1.planning.final)
summary(attr.mod)

plan.nonSC <- subset(lc1.planning.final, Reason.Attendance.Data.Missing!="CCR could not find data")
attr.mod <- lm(attr ~ treat, data=plan.nonSC)
summary(attr.mod)

table(lc1.planning.final$Reason.Attendance.Data.Missing)

par(default.par)
par(mar=c(2.1,4.1,1.1,1.1))
boxplot(Attendance ~ treat, data= lc1.planning, names=c("Control","Treated"),range=0, ylab="Number of Residents Attending Planning Meeting")
means <- tapply(lc1.planning$Attendance,lc1.planning$treat, function(x) mean(x,na.rm=T))
points(means,col="red",pch=18,cex=2)


######################################################################################
###Figure 7: Reanalyzed treatment effects for subgroup of subjects on 
###          neutral and negative side of item scale during the baseline survey
######################################################################################

setwd("~/Google Drive/Bwindi project/Code/WD_Replication")
source('outcome_var_setup.R')
M.ri <- merge(M, sbj_RI, by="Subject_ID_number", all.x=TRUE)

###Model-assisted RI

#Satisfaction with RS Info from UWA
M2.ri <- subset(M.ri, b.knowledge.b7<=3)
mod.h7.ri <- lm.ri(e.knowledge.b7 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                   dta=M2.ri, treat.var = "treat", sims=10000)

#Know how RS works generally
M2.ri <- subset(M.ri, b.knowledge.b8<=2)
mod.h7.2.ri <- lm.ri(e.knowledge.b8 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Know how RS works in my village
M2.ri <- subset(M.ri, b.knowledge.b9<=2)
mod.h7.3.ri <- lm.ri(e.knowledge.b9 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Opportunities to participate in planning
M2.ri <- subset(M.ri, b.participate.b10<=2)
mod.h8.1.ri <- lm.ri(e.participate.b10 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Opportunities to communicate with UWA about RS
M2.ri <- subset(M.ri, b.participate.b11<=3)
mod.h9.1.ri <- lm.ri(e.participate.b11 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Know Right Person to contact about RS
M2.ri <- subset(M.ri, b.participate.b12<=2)
mod.h9.2.ri <- lm.ri(e.participate.b12 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Satisfied with Park Management
M2.ri <- subset(M.ri, b.satisfaction.b2<=3)
nod.h10mgmt.ri <- lm.ri(e.satisfaction.b2 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                        dta=M2.ri, treat.var = "treat", sims=10000)

#Satisfied with RS
M2.ri <- subset(M.ri, b.satisfaction.b3<=3)
nod.h11.ri <- lm.ri(e.satisfaction.b3 ~ treat + b.knowledge.b7 + b.knowledge.b8 + b.knowledge.b9 + b.participate.b10 + b.participate.b11 + b.participate.b12 + b.satisfaction.b2 + b.satisfaction.b3 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                    dta=M2.ri, treat.var = "treat", sims=10000)

#Importance of Protecting Bwindi
#Note: this doesn't work b/c too few observations
#M2.ri <- subset(M.ri, b.conservation.b6<=2)
#nod.h13.ri <- lm.ri(e.conservation.b6 ~ treat + b.conservation.b6 + E_subject_s_gender + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
#                      dta=M2.ri, treat.var = "treat", sims=10000)

#Plot
names <- c("Satisfied with Park Management","Satisfied with Revenue Sharing","Satisfied with Information from UWA about RS","Opportunities to Participate in RS","Opportunities to Communicate with UWA about RS","Know Right Person to Contact about RS","Know How RS Works in my Village","Know How RS Works Generally")
ates <- c(nod.h10mgmt.ri$ate,nod.h11.ri$ate,mod.h7.ri$ate,mod.h8.1.ri$ate,mod.h9.1.ri$ate,mod.h9.2.ri$ate,mod.h7.3.ri$ate,mod.h7.2.ri$ate)
ses <- c(nod.h10mgmt.ri$se,nod.h11.ri$se,mod.h7.ri$se,mod.h8.1.ri$se,mod.h9.1.ri$se,mod.h9.2.ri$se,mod.h7.3.ri$se,mod.h7.2.ri$se)
dta.coef <- data.frame(names,ates,ses)
dta.coef$ylo <- dta.coef$ates - 1.645*dta.coef$ses
dta.coef$yhi <- dta.coef$ates + 1.645*dta.coef$ses
dta.coef$names <- factor(dta.coef$names, levels=dta.coef$names)
dta.coef$clr <- c(rep("blue",3),rep("red",2),rep("darkgreen",3))

p2.low = ggplot(dta.coef, aes(x=names, y=ates, ymin=ylo, ymax=yhi)) + geom_errorbar(colour=dta.coef$clr, width=0.2, size=1.4) + 
  theme_bw(base_size=16) + coord_flip() + geom_hline(aes(yintercept=0), lty=2) + ylab("Treatment Effects on Survey Items") + xlab("") + geom_point(size=4, colour=dta.coef$clr) +
  theme(axis.title = element_text(size=14))
p2.low #Copied at 640x420


######################################################################################
###Figure 8: Reanalyzed treatment effects for women 
######################################################################################

setwd("~/Google Drive/Bwindi project/Code/WD_Replication")
source('outcome_var_mean_imp_setup.R')
M.ri <- merge(M, sbj_RI, by="Subject_ID_number", all.x=TRUE)

M2.ri <- subset(M.ri, E_subject_s_gender=="female")

#Satisfaction with RS Info from UWA
mod.h7.ri <- lm.ri(e.knowledge.b7 ~ treat + b.knowledge.b7 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                   dta=M2.ri, treat.var = "treat", sims=10000)

#Know how RS works generally
mod.h7.2.ri <- lm.ri(e.knowledge.b8 ~ treat + b.knowledge.b8 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Know how RS works in my village
mod.h7.3.ri <- lm.ri(e.knowledge.b9 ~ treat + b.knowledge.b9 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Opportunities to participate in planning
mod.h8.1.ri <- lm.ri(e.participate.b10 ~ treat + b.participate.b10 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Participate in RS meetings in past months
mod.h8.3.ri <- lm.ri(e.participate.e4 ~ treat + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)
#Note: this result does not match exactly the result reported in Figure 8 due to an error uncovered during the preparation of replication files

#Opportunities to communicate with UWA about RS
mod.h9.1.ri <- lm.ri(e.participate.b11 ~ treat + b.participate.b11 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Know Right Person to contact about RS
mod.h9.2.ri <- lm.ri(e.participate.b12 ~ treat + b.participate.b12 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                     dta=M2.ri, treat.var = "treat", sims=10000)

#Satisfied with Park Management
nod.h10mgmt.ri <- lm.ri(e.satisfaction.b2 ~ treat + b.satisfaction.b2 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                        dta=M2.ri, treat.var = "treat", sims=10000)

#Satisfied with RS
nod.h11.ri <- lm.ri(e.satisfaction.b3 ~ treat + b.satisfaction.b3 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                    dta=M2.ri, treat.var = "treat", sims=10000)

#Importance of Protecting Bwindi
nod.h13.ri <- lm.ri(e.conservation.b6 ~ treat + b.conservation.b6 + E_Subject_s_age + E_Subject_s_income + E_Can_you_read_001 + subjects.per.village + subcounty.blocking,
                    dta=M2.ri, treat.var = "treat", sims=10000)

#Plot
names <- c("Importance of Protecting Bwindi","Satisfied with Park Management","Satisfied with Revenue Sharing","Satisfied with Information from UWA about RS","Participated in RS meeting (past month)","Opportunities to Participate in RS","Opportunities to Communicate with UWA about RS","Know Right Person to Contact about RS","Know How RS Works in my Village","Know How RS Works Generally")
ates <- c(nod.h13.ri$ate,nod.h10mgmt.ri$ate,nod.h11.ri$ate,mod.h7.ri$ate,mod.h8.3.ri$ate,mod.h8.1.ri$ate,mod.h9.1.ri$ate,mod.h9.2.ri$ate,mod.h7.3.ri$ate,mod.h7.2.ri$ate)
ses <- c(nod.h13.ri$se,nod.h10mgmt.ri$se,nod.h11.ri$se,mod.h7.ri$se,mod.h8.3.ri$se,mod.h8.1.ri$se,mod.h9.1.ri$se,mod.h9.2.ri$se,mod.h7.3.ri$se,mod.h7.2.ri$se)
dta.coef <- data.frame(names,ates,ses)
dta.coef$ylo <- dta.coef$ates - 1.645*dta.coef$ses
dta.coef$yhi <- dta.coef$ates + 1.645*dta.coef$ses
dta.coef$names <- factor(dta.coef$names, levels=dta.coef$names)
dta.coef$clr <- c(rep("blue",4),rep("red",3),rep("darkgreen",3))

p2.women = ggplot(dta.coef, aes(x=names, y=ates, ymin=ylo, ymax=yhi)) + geom_errorbar(colour=dta.coef$clr, width=0.2, size=1.4) + 
  theme_bw(base_size=16) + coord_flip() + geom_hline(aes(yintercept=0), lty=2) + ylab("Treatment Effects on Survey Items") + xlab("") + geom_point(size=4, colour=dta.coef$clr) +
  theme(axis.title = element_text(size=14))
p2.women #Copied at 800x420


######################################################################################
###Figure 9: Conditional treatment effects by number of subjects reached 
######################################################################################

setwd("~/Google Drive/Bwindi project/Code/WD_Replication")
source('outcome_var_mean_imp_setup.R')
M.ri <- merge(M, sbj_RI, by="Subject_ID_number", all.x=TRUE)

sat.bin <- inter.binning.90(Y="e.satisfaction.b2", D="treat", X="subjects.per.village", Z="b.satisfaction.b2", FE="subcounty.blocking", data=M, 
                            Ylabel = "Outcome", Dlabel = "Treatment", Xlabel="Number of Subjects Per Village", main="Satisfied with Park Management", nbins=4, na.rm=TRUE,
                            ylim=c(-0.72,0.85))
park.sat.bin <- sat.bin$graph

info.bin <- inter.binning.90(Y="e.knowledge.b7", D="treat", X="subjects.per.village", Z="b.knowledge.b7", FE="subcounty.blocking", data=M, 
                             Ylabel = "Outcome", Dlabel = "Treatment", Xlabel="Number of Subjects Per Village", main="Satisfied with Information from UWA about RS", nbins=4, na.rm=TRUE)
info.from.UWA.bin <- info.bin$graph

import.bin <- inter.binning.90(Y="e.conservation.b6", D="treat", X="subjects.per.village", Z="b.conservation.b6", FE="subcounty.blocking", data=M, 
                               Ylabel = "Outcome", Dlabel = "Treatment", Xlabel="Number of Subjects Per Village", main="Importance of Protecting Bwindi", nbins=4, na.rm=TRUE)
con.support.bin <- import.bin$graph

partic.bin <- inter.binning.90(Y="e.participate.e4", D="treat", X="subjects.per.village", FE="subcounty.blocking", data=M, 
                               Ylabel = "Outcome", Dlabel = "Treatment", Xlabel="Number of Subjects Per Village", main="Participated in RS Meeting (past month)", nbins=4, na.rm=TRUE)
part.bin <- partic.bin$graph

#pdf("BinningMEs_170928.pdf", width=8.5, height=6)
multiplot(part.bin,park.sat.bin,info.from.UWA.bin,con.support.bin, cols=2) #For binning estimates, WD.R1
#dev.off()



